#Table 4 (in Suppl materials): Alternative dynamic spatial autocorrelation (SAC) models

#Model 1 (1975-2021)
#####################################################################################

# we need original data 'medata_1975_2021' 
data <- medata_1975_2021
data = data[order(data$year, data$ccode),]
usa = subset(data, data$ccode==2)
data=subset(data, data$ccode !=2) #data without USA

#colnames(usa)
names(usa)[names(usa) == "me_abs"] <- "usa.me_abs"
names(usa)[names(usa) == "me_abslag"] <- "usa.me_abslag"
names(usa)[names(usa) == "lnme_abslag"] <- "usa.lnme_abslag"
colnames(usa)
usa <- subset(usa, select = c("country", "year", "usa.me_abs", "usa.me_abslag", "usa.lnme_abslag"))
#colnames(usa)
#[1] "country"         "year"            "usa.me_abs"      "usa.me_abslag"   "usa.lnme_abslag"

data1=merge(data, usa, by="year") 

data.ordered <- data1[order(data1$year, data1$ccode),] #to be sure of appropriate ordering
#View(data.ordered)

#download corresponding spatial matrix without USA
mematrix1975_2021_without_usa <- read_excel("mematrix1975_2021_without_usa.xlsx")
capdistinv.wide=as.matrix(mematrix1975_2021_without_usa)

capdistinv.wide.resc <- capdistinv.wide*100 #re-scale inverse distances in hundreds of km
capdistinv.2listw.nors.resc <- mat2listw(capdistinv.wide.resc) #non-normalized (re-scaled) weights matrix
##'Warning message' because non-normalized

listw.1975 = capdistinv.2listw.nors.resc

#counters
counter=seq(1,47,1)
y=seq(1975,2021,1)
counter=as.matrix(cbind(y, counter))
colnames(counter) = c("year","counter")
data.ordered=merge(data.ordered, counter, by="year")
data.ordered$counter.sq=data.ordered$counter^2
data.ordered$counter.cub=data.ordered$counter^3

#without US
sacL.resc.no.us <- sacsarlm(lnme_abs ~ lnme_abslag  + lngdp_abs + lnpop + as.factor(ccode)+counter+
                              counter.sq+counter.cub, data = data.ordered, listw.1975, Durbin = F)
sacL.resc.no.us$rho
#summary(sacL.resc.no.us)
#####################################################################################


#Model 2 (1993-2021 & rusp_(t-1))
#####################################################################################

data1993 = subset(medata_1975_2021, medata_1975_2021$year > 1992)
data1993=subset(data1993, data1993$ccode !=2) #data without USA
data.ordered1993 = data1993[order(data1993$year, data1993$ccode),] 
dim(data.ordered1993)

#download corresponding matrix (1993-2021)
mematrix1993_2021_rusp_dist_without_usa <- read_excel("mematrix1993_2021_rusp_dist_without usa.xlsx")
capdistinv.wide <- as.matrix(mematrix1993_2021_rusp_dist_without_usa)
dim(capdistinv.wide) 

capdistinv.wide.resc <- capdistinv.wide*100 #re-scale inverse distances in terms of hundreds of km
capdistinv.2lw.nors.resc <- mat2listw(capdistinv.wide.resc) #non-normalized (re-scaled) weights matrix

listw.1993.ru = capdistinv.2lw.nors.resc

#counter (time trends)
length(unique(data.ordered1993$year))
counter=seq(1,29,1)
y=seq(1993,2021,1)
counter=as.matrix(cbind(y, counter))
colnames(counter) = c("year","counter")
data.ordered1993=merge(data.ordered1993, counter, by="year")
data.ordered1993$counter.sq=data.ordered1993$counter^2
data.ordered1993$counter.cub=data.ordered1993$counter^3

data.ordered1993 = data.ordered1993[order(data.ordered1993$year, data.ordered1993$ccode),] #be sure of appropriate ordering for spatial analysis

#Model 2 (1993-2021), with russp (lagged Russian spending, ln)
sacL.resc1993.ru <- sacsarlm(lnme_abs ~ lnme_abslag + lngdp_abs + lnpop + lLru_abs + as.factor(ccode) + counter+counter.sq + 
                               counter.cub, data = data.ordered1993, listw.1993.ru, Durbin = F) #weights non-normalized
sacL.resc1993.ru$rho
#summary(sacL.resc1993.ru)
#####################################################################################


#Model 3 (1993-2021), with russp (lagged, ln) and distance to RU
#####################################################################################

sacL.resc1993.distru <- sacsarlm(lnme_abs ~ lnme_abslag + lngdp_abs + lnpop + lLru_abs + lLru_abs*dist_ru+ as.factor(ccode) + 
                                   counter+counter.sq + counter.cub, data = data.ordered1993, listw.1993.ru, Durbin = F) #weights non-normalized
                                   ## 'Warning message' because of 'dust_ru' presence
sacL.resc1993.distru$rho
#####################################################################################

###Output
###Models 1-3

#install.package("texreg")
#library(texreg)
#capture.output(screenreg(list(sacL.resc.no.us, sacL.resc1993.ru, sacL.resc1993.distru), stars = c(0.1, 0.05, 0.01), digits=3), file = "output_Table4_Models1-3_in_SM.txt")


#Model 4 (with US as nearest neighbor)
#####################################################################################

data2 = subset(medata_1975_2021, medata_1975_2021$year > 1992) 
data.ordered = data2[order(data2$year, data2$ccode),]
dim(data.ordered) #652 27

#download corresponding matrix
mematrix1993_2021_US_as_nearest_neigh <- read_excel("mematrix1993-2021_US as nearest neigh.xlsx")
dist.wide = as.matrix(mematrix1993_2021_US_as_nearest_neigh)
dim(dist.wide) #[1] 652 652

dist.lw <- mat2listw(dist.wide) #non-normalized (re-scaled) weights matrix
lw.dist.1993 = dist.lw

colnames(data.ordered)
#counter
length(unique(data.ordered$year))
counter=seq(1,29,1)
y=seq(1993,2021,1)
counter=as.matrix(cbind(y, counter))
colnames(counter) = c("year","counter")
data.ordered=merge(data.ordered, counter, by="year")
data.ordered$counter.sq=data.ordered$counter^2
data.ordered$counter.cub=data.ordered$counter^3
dim(data.ordered)

data.ordered = data.ordered[order(data.ordered$year, data.ordered$ccode),] #be sure of appropriate ordering for spatial analysis

sacL.dist.usneigh <- sacsarlm(lnme_abs ~ lnme_abslag  + lngdp_abs + lnpop + lLru_abs + lLru_abs * dist_ru + as.factor(ccode)+
                                counter+counter.sq+ counter.cub,data = data.ordered, dist.lw, Durbin = F, zero.policy = TRUE)
sacL.dist.usneigh$rho
#summary(sacL.dist.usneigh)

#capture.output(screenreg(list(sacL.dist.usneigh), stars = c(0.1, 0.05, 0.01), digits=3), file = "Table4_Model4_in_SM.txt")
#####################################################################################


#Model 5 (Host-nation support)
#####################################################################################

#with HNS index 
#source: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/ESQQAG ('ReplicationData')
#when merged with ME data, we have 'ReplicationData_merged'
ReplicationData_merged <- read_excel("ReplicationData_ merged.xlsx")
rep=ReplicationData_merged

#sort(unique(rep$year))
#[1] 1995 1996 1997 1998 1999 2000 2001 2002

rep$country[rep$country == "Slovak Republic"] = "Slovakia"

rep = rep[order(rep$year, rep$ccode),]
rep$year_country <- rep$year*1000 + rep$ccode
data.dfi <- rep%>% select("year_country") ## one needs dplyr for this
data.dfj <- data.dfi
colnames(data.dfi) <- "year_countryi"
colnames(data.dfj) <- "year_countryj"
data_cross <- merge(data.dfi, data.dfj, all=TRUE)
#View(data_cross)

#to attribute a ccodeij identifier for each country pair (will need for attributing values on distances)
data_cross$yeari <- floor(data_cross$year_countryi/1000)
data_cross$yearj <- floor(data_cross$year_countryj/1000)
data_cross$ccodei <- data_cross$year_countryi - 1000*data_cross$yeari
data_cross$ccodej <- data_cross$year_countryj - 1000*data_cross$yearj
data_cross$ccodeij <- data_cross$ccodei*1000 + data_cross$ccodej 
# View(data_cross) 

# we have "long" data.frame (another file) with distances and ccodeij identifier for each country pair 
long <- read_excel("long.xlsx")

# loop to fill data - attributes values on distances to each country pair (across years)
for (i in 1:nrow(data_cross)) {
  data_cross$capdist[i] <- long$capdist[long$ccodeij == data_cross$ccodeij[i]]
}
#View(data_cross)

# no association across time 
capdist_cross <- data_cross
capdist_cross$capdist <- ifelse(capdist_cross$yeari != capdist_cross$yearj, 0, capdist_cross$capdist)
#View(capdist_cross)

# inverse
capdist_cross$capdist_inv <- ifelse(capdist_cross$capdist == 0, 0, 1/capdist_cross$capdist)
#View(capdist_cross)
capdist <- subset(capdist_cross, select=c("year_countryi", "year_countryj", "capdist_inv"))

## Connectivity matrix
capdistinv.wide <- reshape(capdist, timevar="year_countryj", idvar="year_countryi", direction="wide") 
#View(capdistinv.wide)
#dim(capdistinv.wide)
capdistinv.wide <- as.data.frame(capdistinv.wide[, -1]) ## drop id column
dim(capdistinv.wide) 

capdistinv.wide=as.matrix(capdistinv.wide)
capdistinv.wide.resc <- capdistinv.wide*100 #re-scaled
capdistinv.2listw.nors.resc <- mat2listw(capdistinv.wide.resc) #non-normalized re-scaled weights matrix

listw.1995.rep = capdistinv.2listw.nors.resc

#counter
length(unique(rep$year))
counter=seq(1,8,1)
y=seq(1995,2002,1)
counter=as.matrix(cbind(y, counter))
colnames(counter) = c("year","counter")
rep=merge(rep, counter, by="year")
rep$counter.sq=rep$counter^2
rep$counter.cub=rep$counter^3

#sac
sacL.rep1995 <- sacsarlm(lnme_abs ~ lnme_abslag + lngdp_abs + lnpop + hns_pctgdp3 + lLru_abs + lLru_abs * dist_ru +
                           as.factor(ccode) + counter+counter.sq+counter.cub, data = rep, listw.1995.rep, Durbin = F) #weights non-normalized

#'Warning message' regarding the use of an alternative method for the inversion of the asymptotic covariance matrix may be due to a small sample size     
##dropping covariate 'distr_ru' may serve as a robustness check for the results of interest - the spatial lag and HNS coefficients 

sacL.rep1995$rho
#summary(sacL.rep1995)

#capture.output(screenreg(list(sacL.rep1995), stars = c(0.1, 0.05, 0.01), digits=3), file = "Table4_Model5_in_SM.txt")

capture.output(screenreg(list(sacL.resc.no.us, sacL.resc1993.ru, sacL.resc1993.distru,sacL.dist.usneigh,sacL.rep1995), 
                         stars = c(0.1, 0.05, 0.01), digits=3), file = "Table4_SM.txt")

#####################################################################################

